Na naszym ostatnim spotkaniu Tomek pokazywał wykresy z falek z których wynikało, że jest taki przedział CSI w którym wydaje się że jest sprzężenie (i to istotne statystycznie).
# load data
low_data_dir <- "../../data/low_resolution"
high_data_dir <- "../../data/high_resolution"
low_files <- dir(low_data_dir, patter = "*.csv$")
high_files <- dir(high_data_dir, patter = "*.csv$")
low_data <- low_files %>%
map(~ read_csv(file.path(low_data_dir, .), show_col_types = F)) %>%
reduce(rbind)
high_data <- high_files %>%
map(~ read_csv(file.path(high_data_dir, .), show_col_types = F)) |>
reduce(rbind)
# high res data gathered in 4 sessions encoded PARTID + S1-S4. S1-S4 removed
high_data <- high_data |>
mutate(PART_ID = str_sub(PART_ID, 1, nchar(PART_ID)-2))
Dane surowe (high freq)
high_data %>%
group_by(PART_ID, CSI) %>%
summarise(mean = mean(Corr)) %>%
ggplot(mapping = aes(x = CSI, y = mean)) +
geom_line() +
facet_wrap(~PART_ID)

Teoria
Transformata fouriera działa tak, że bierzemy sobie sygnał i rozbijamy go na sumę cosinusów o różnej częstotliwości i fazie. Wiemy też, że transformata fouriera działa w obie strony - na zbiorze cosinusów możemy zrobić transformatę odwrotną inverse-FFT i dostajemy sygnał oryginalny. Z kolei z falkami jest tak, że zamiast cosinusów mamy jedną funkcję bazową, najpopularniejszym wyborem jest falka morleta; 
(Formalnie to trzeba by było pamiętać, że tam jest komponent rzeczywisty i urojony, ale dla wywodu to nie ma znaczenia)
No i transformata falkowa to jest takie zwierze, że bierze tyle takich falek morleta ile potrzebuje i każdą z nich rozciąga/ściska a następnie przesuwa po osi x. Czyli tak jak wynikiem Fouriera jest zbiór współczynników jak przeskalować (jaką amplitudę nadać) i jak przesunąć w fazie (ale o tym na potrzebę wywodu zapominamy, bo to część urojona) kolejne cosinusy, to wynikiem falek jest recepta jak pougniatać i przeciągać po osi x ileś tam kopii funkcji bazowej.
No i tutaj się pojawia sanity-check w postaci rekonstrukcji sygnału. Jest to nic innego jak zrobienie odwrotnej transformaty falkowej na tym zbiorze pogniecionych i przesuwanych falek w celu uzyskania oryginalnego sygnału.
I to właśnie widać na wykresach poniżej.
Wykresy
Przepraszam, że zostawiam brzydko opisane ośki. Biblioteka która rysuje te falki średnio ze mną współpracuje i uznałem, że użeranie się z nią nie ma sensu.
Dla każdego badanego mamy trzy wykresy, kolejno:
- Widmo czasowo-częstościowe
- oś X to kolejne punkty danych, czyli ten nasz przedział 300-900 ms czy jakoś tak.
- oś Y - długość okresu czyli 1/Hz skala logarytmiczna 0.25 to 4 Hz 0.125 to 8 Hz 0.0625 to 16 Hz 0.03125 to 32 Hz
- Białą obwolutą zaznaczone są obszary istotne statystycznie
- Falki
Poszczególne falki dopasowane do sygnału. Warte uwagi jest to, że (1) jest ich relatywnie niewiele (2) są duże obszary pustki, co oznacza, że aproksymujemy zerem.
- Wynik rekonstrukcji falek nałożony na oryginalny sygnał.
Widać, że te rekonstrukcje są kiepskie, co oznacza, że dopasowanie falek jest równie złe.
Insight: Właściwie wszystkie wykresy mają obszary istotne statystycznie, zgodnie z wynikami Tomka. Są one jednak rozmieszczone dość chaotycznie.
Insight: Ta istotność statystyczna wychodzi być może z porównywania obszarów gdzie amplituda dopasowanych falek jest w miare duża, do obszarów, w których nie dopasowano niczego. Z porównania wykresów widm do rekonstrukcji widać, że istotnie jest tam, gdzie jest peak falki.
Insight: Te “wyspy istotności statystycznej” są bardzo niestabilne numerycznie. Wystarczy lekko zmienić skale osi X (nawet na wielokrotność) i wynik jest zupełnie inny. Gdyby te efekty były realne, to zmiana parametru próbkowania nie powinna wiele zmieniać, co najwyżej zakłamywać odczyt lokalizacji z wykresu, tutaj natomiast po zmianie parametru, “wyspy” pojawiają się zupełnie gdzie indziej. To również sugeruje, że istotność może wychodzić z numerycznych błędów przy porównywaniu do zera.
part_ids <- high_data |> select(PART_ID) |> unique()
part_ids <- part_ids$PART_ID
# params
detrend = FALSE
sampling = 120
low_freq_band = 4
upper_freq_band = 32
for (pid in part_ids) {
w <- high_data |>
filter(Trial_type == "experiment", PART_ID == pid) |>
group_by(CSI) |>
summarise(mean = mean(Corr)) |>
analyze.wavelet(
"mean",
loess.span = detrend,
dt = 1 / sampling,
lowerPeriod = 1 / upper_freq_band,
upperPeriod = 1 / low_freq_band,
make.pval = TRUE,
n.sim = 10
)
wt.image(
w,
color.key = "quantile",
main = pid,
n.levels = 100,
legend.params = list(lab = "wavelet power levels", mar = 4.7)
)
reconstruct(
w,
plot.waves = TRUE,
lwd = c(1, 2),
legend.coords = "bottomleft",
main = pid
)
}



























LS0tCnRpdGxlOiAiYHIgcGFyYW1zJGR5bmFtaWN0aXRsZWAiCmF1dGhvcjogIkJhcnTFgm9taWVqIEtyb2N6ZWsiCmRhdGU6ICJgciBTeXMuRGF0ZSgpYCIKcGFyYW1zOgogIGR5bmFtaWN0aXRsZTogV2F2ZWxldHMgZm9yIGJlaGF2aW9yYWwgb3NjaWxsYXRpb25zCiAgdmlyaWRpc19wYWxldHRlOiB2aXJpZGlzCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdG9jOiB0cnVlCiAgICB0b2NfZmxvYXQ6IGZhbHNlCiAgICBudW1iZXJfc2VjdGlvbnM6IGZhbHNlCiAgICB0aGVtZTogImZsYXRseSIKICAgIGhpZ2hsaWdodDogdGV4dG1hdGUKICAgIGNvZGVfZm9sZGluZzogc2hvdwogICAgY29kZV9kb3dubG9hZDogVFJVRQotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlID0gRkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldCh3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSkKZ2dwbG90Mjo6dGhlbWVfc2V0KGdncGxvdDI6OnRoZW1lX21pbmltYWwoKSkKYGBgCgpgYGB7ciBpbmNsdWRlPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShXYXZlbGV0Q29tcCkKYGBgCgpOYSBuYXN6eW0gb3N0YXRuaW0gc3BvdGthbml1IFRvbWVrIHBva2F6eXdhxYIgd3lrcmVzeSB6IGZhbGVrIHoga3TDs3J5Y2ggd3luaWthxYJvLCDFvGUgamVzdCB0YWtpIHByemVkemlhxYIgQ1NJIHcga3TDs3J5bSB3eWRhamUgc2nEmSDFvGUgamVzdCBzcHJ6xJnFvGVuaWUgKGkgdG8gaXN0b3RuZSBzdGF0eXN0eWN6bmllKS4KCmBgYHtyIHJlc3VsdHM9RkFMU0V9CiMgbG9hZCBkYXRhCmxvd19kYXRhX2RpciA8LSAiLi4vLi4vZGF0YS9sb3dfcmVzb2x1dGlvbiIKaGlnaF9kYXRhX2RpciA8LSAiLi4vLi4vZGF0YS9oaWdoX3Jlc29sdXRpb24iCmxvd19maWxlcyA8LSBkaXIobG93X2RhdGFfZGlyLCBwYXR0ZXIgPSAiKi5jc3YkIikKaGlnaF9maWxlcyA8LSBkaXIoaGlnaF9kYXRhX2RpciwgcGF0dGVyID0gIiouY3N2JCIpCgpsb3dfZGF0YSA8LSBsb3dfZmlsZXMgJT4lCiAgbWFwKH4gcmVhZF9jc3YoZmlsZS5wYXRoKGxvd19kYXRhX2RpciwgLiksIHNob3dfY29sX3R5cGVzID0gRikpICU+JQogIHJlZHVjZShyYmluZCkKCmhpZ2hfZGF0YSA8LSBoaWdoX2ZpbGVzICU+JQogIG1hcCh+IHJlYWRfY3N2KGZpbGUucGF0aChoaWdoX2RhdGFfZGlyLCAuKSwgc2hvd19jb2xfdHlwZXMgPSBGKSkgfD4KICByZWR1Y2UocmJpbmQpCgojIGhpZ2ggcmVzIGRhdGEgZ2F0aGVyZWQgaW4gNCBzZXNzaW9ucyBlbmNvZGVkIFBBUlRJRCArIFMxLVM0LiBTMS1TNCByZW1vdmVkCmhpZ2hfZGF0YSA8LSBoaWdoX2RhdGEgfD4gCiAgICAgICAgICAgICAgbXV0YXRlKFBBUlRfSUQgPSBzdHJfc3ViKFBBUlRfSUQsIDEsIG5jaGFyKFBBUlRfSUQpLTIpKQpgYGAKCgpgYGB7ciBldmFsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQpzdHIoaGlnaF9kYXRhKQpgYGAKCiMjIERhbmUgc3Vyb3dlIChoaWdoIGZyZXEpCgpgYGB7ciBmaWcud2lkdGg9MTJ9CmhpZ2hfZGF0YSAlPiUKICBncm91cF9ieShQQVJUX0lELCBDU0kpICU+JQogIHN1bW1hcmlzZShtZWFuID0gbWVhbihDb3JyKSkgJT4lCiAgZ2dwbG90KG1hcHBpbmcgPSBhZXMoeCA9IENTSSwgeSA9IG1lYW4pKSArCiAgZ2VvbV9saW5lKCkgKwogIGZhY2V0X3dyYXAoflBBUlRfSUQpCmBgYAoKIyMgVGVvcmlhIAoKVHJhbnNmb3JtYXRhIGZvdXJpZXJhIGR6aWHFgmEgdGFrLCDFvGUgYmllcnplbXkgc29iaWUgc3lnbmHFgiBpIHJvemJpamFteSBnbyBuYSBzdW3EmSBjb3NpbnVzw7N3IG8gcsOzxbxuZWogY3rEmXN0b3RsaXdvxZtjaSBpIGZhemllLgpXaWVteSB0ZcW8LCDFvGUgdHJhbnNmb3JtYXRhIGZvdXJpZXJhIGR6aWHFgmEgdyBvYmllIHN0cm9ueSAtIG5hIHpiaW9yemUgY29zaW51c8OzdyBtb8W8ZW15IHpyb2JpxIcgdHJhbnNmb3JtYXTEmSBvZHdyb3RuxIUgaW52ZXJzZS1GRlQgaSBkb3N0YWplbXkgc3lnbmHFgiBvcnlnaW5hbG55LgpaIGtvbGVpIHogZmFsa2FtaSBqZXN0IHRhaywgxbxlIHphbWlhc3QgY29zaW51c8OzdyBtYW15IGplZG7EhSBmdW5rY2rEmSBiYXpvd8SFLCBuYWpwb3B1bGFybmllanN6eW0gd3lib3JlbSBqZXN0IGZhbGthIG1vcmxldGE7CiFbTW9ybGV0IHdhdmVsZXQgc2hhcGVdKGltZy9tb3JsZXQucG5nKXt3aWR0aD01MCV9CgooRm9ybWFsbmllIHRvIHRyemViYSBieSBiecWCbyBwYW1pxJl0YcSHLCDFvGUgdGFtIGplc3Qga29tcG9uZW50IHJ6ZWN6eXdpc3R5IGkgdXJvam9ueSwgYWxlIGRsYSB3eXdvZHUgdG8gbmllIG1hIHpuYWN6ZW5pYSkKCgpObyBpIHRyYW5zZm9ybWF0YSBmYWxrb3dhIHRvIGplc3QgdGFraWUgendpZXJ6ZSwgxbxlIGJpZXJ6ZSB0eWxlIHRha2ljaCBmYWxlayBtb3JsZXRhIGlsZSBwb3RyemVidWplIGkga2HFvGTEhSB6IG5pY2ggcm96Y2nEhWdhL8WbY2lza2EgYSBuYXN0xJlwbmllIHByemVzdXdhIHBvIG9zaSB4LiBDenlsaSB0YWsgamFrIHd5bmlraWVtIEZvdXJpZXJhIGplc3QgemJpw7NyIHdzcMOzxYJjenlubmlrw7N3IGphayBwcnplc2thbG93YcSHIChqYWvEhSBhbXBsaXR1ZMSZIG5hZGHEhykgaSBqYWsgcHJ6ZXN1bsSFxIcgdyBmYXppZSAoYWxlIG8gdHltIG5hIHBvdHJ6ZWLEmSB3eXdvZHUgemFwb21pbmFteSwgYm8gdG8gY3rEmcWbxIcgdXJvam9uYSkga29sZWpuZSBjb3NpbnVzeSwgdG8gd3luaWtpZW0gZmFsZWsgamVzdCByZWNlcHRhIGphayBwb3VnbmlhdGHEhyBpIHByemVjacSFZ2HEhyBwbyBvc2kgeCBpbGXFmyB0YW0ga29waWkgZnVua2NqaSBiYXpvd2VqLgoKTm8gaSB0dXRhaiBzacSZIHBvamF3aWEgc2FuaXR5LWNoZWNrIHcgcG9zdGFjaSByZWtvbnN0cnVrY2ppIHN5Z25hxYJ1LiBKZXN0IHRvIG5pYyBpbm5lZ28gamFrIHpyb2JpZW5pZSBvZHdyb3RuZWogdHJhbnNmb3JtYXR5IGZhbGtvd2VqIG5hIHR5bSB6YmlvcnplIHBvZ25pZWNpb255Y2ggaSBwcnplc3V3YW55Y2ggZmFsZWsgdyBjZWx1IHV6eXNrYW5pYSBvcnlnaW5hbG5lZ28gc3lnbmHFgnUuIAoKSSB0byB3xYJhxZtuaWUgd2lkYcSHIG5hIHd5a3Jlc2FjaCBwb25pxbxlai4KCiMjIFd5a3Jlc3kKClByemVwcmFzemFtLCDFvGUgem9zdGF3aWFtIGJyenlka28gb3Bpc2FuZSBvxZtraS4gQmlibGlvdGVrYSBrdMOzcmEgcnlzdWplIHRlIGZhbGtpIMWbcmVkbmlvIHplIG1uxIUgd3Nww7PFgnByYWN1amUgCmkgdXpuYcWCZW0sIMW8ZSB1xbxlcmFuaWUgc2nEmSB6IG5pxIUgbmllIG1hIHNlbnN1LiAKCkRsYSBrYcW8ZGVnbyBiYWRhbmVnbyBtYW15IHRyenkgd3lrcmVzeSwga29sZWpubzoKCjEuIFdpZG1vIGN6YXNvd28tY3rEmXN0b8WbY2lvd2UKCiogb8WbIFggdG8ga29sZWpuZSBwdW5rdHkgZGFueWNoLCBjenlsaSB0ZW4gbmFzeiBwcnplZHppYcWCIDMwMC05MDAgbXMgY3p5IGpha2/FmyB0YWsuIAoqIG/FmyBZIC0gZMWCdWdvxZvEhyBva3Jlc3UgY3p5bGkgMS9IeiBza2FsYSBsb2dhcnl0bWljem5hIDAuMjUgdG8gNCBIeiAwLjEyNSB0byA4IEh6IDAuMDYyNSB0byAxNiBIeiAwLjAzMTI1IHRvIDMyIEh6CiogQmlhxYLEhSBvYndvbHV0xIUgemF6bmFjem9uZSBzxIUgb2JzemFyeSBpc3RvdG5lIHN0YXR5c3R5Y3puaWUKCjIuIEZhbGtpCgpQb3N6Y3plZ8OzbG5lIGZhbGtpIGRvcGFzb3dhbmUgZG8gc3lnbmHFgnUuIFdhcnRlIHV3YWdpIGplc3QgdG8sIMW8ZSAoMSkgamVzdCBpY2ggcmVsYXR5d25pZSBuaWV3aWVsZSAoMikgc8SFIGR1xbxlIG9ic3phcnkKcHVzdGtpLCBjbyBvem5hY3phLCDFvGUgYXByb2tzeW11amVteSB6ZXJlbS4gCgozLiBXeW5payByZWtvbnN0cnVrY2ppIGZhbGVrIG5hxYJvxbxvbnkgbmEgb3J5Z2luYWxueSBzeWduYcWCLiAKCldpZGHEhywgxbxlIHRlIHJla29uc3RydWtjamUgc8SFIGtpZXBza2llLCBjbyBvem5hY3phLCDFvGUgZG9wYXNvd2FuaWUgZmFsZWsgamVzdCByw7N3bmllIHrFgmUuIAoKKipJbnNpZ2h0Kio6IFfFgmHFm2Npd2llIHdzenlzdGtpZSB3eWtyZXN5IG1hasSFIG9ic3phcnkgaXN0b3RuZSBzdGF0eXN0eWN6bmllLCB6Z29kbmllIHogd3luaWthbWkgVG9ta2EuIFPEhSBvbmUgamVkbmFrIHJvem1pZXN6Y3pvbmUgZG/Fm8SHIGNoYW90eWN6bmllLgoKKipJbnNpZ2h0Kio6IFRhIGlzdG90bm/Fm8SHIHN0YXR5c3R5Y3puYSB3eWNob2R6aSBiecSHIG1vxbxlIHogcG9yw7N3bnl3YW5pYSBvYnN6YXLDs3cgZ2R6aWUgYW1wbGl0dWRhIGRvcGFzb3dhbnljaCBmYWxlayBqZXN0IHcgbWlhcmUgZHXFvGEsIGRvIG9ic3phcsOzdywgdyBrdMOzcnljaCBuaWUgZG9wYXNvd2FubyBuaWN6ZWdvLiBaIHBvcsOzd25hbmlhIHd5a3Jlc8OzdyB3aWRtIGRvIHJla29uc3RydWtjamkgd2lkYcSHLCDFvGUgaXN0b3RuaWUgamVzdCB0YW0sIGdkemllIGplc3QgcGVhayBmYWxraS4gCgoqKkluc2lnaHQqKjogVGUgInd5c3B5IGlzdG90bm/Fm2NpIHN0YXR5c3R5Y3puZWoiIHPEhSBiYXJkem8gbmllc3RhYmlsbmUgbnVtZXJ5Y3puaWUuIFd5c3RhcmN6eSBsZWtrbyB6bWllbmnEhyBza2FsZSBvc2kgWCAobmF3ZXQgbmEgd2llbG9rcm90bm/Fm8SHKSBpIHd5bmlrIGplc3QgenVwZcWCbmllIGlubnkuIEdkeWJ5IHRlIGVmZWt0eSBiecWCeSByZWFsbmUsIHRvIHptaWFuYSBwYXJhbWV0cnUgcHLDs2Jrb3dhbmlhCm5pZSBwb3dpbm5hIHdpZWxlIHptaWVuaWHEhywgY28gbmFqd3nFvGVqIHpha8WCYW15d2HEhyBvZGN6eXQgbG9rYWxpemFjamkgeiB3eWtyZXN1LCB0dXRhaiBuYXRvbWlhc3QgcG8gem1pYW5pZSBwYXJhbWV0cnUsIAoid3lzcHkiIHBvamF3aWFqxIUgc2nEmSB6dXBlxYJuaWUgZ2R6aWUgaW5kemllai4gVG8gcsOzd25pZcW8IHN1Z2VydWplLCDFvGUgaXN0b3Rub8WbxIcgbW/FvGUgd3ljaG9kemnEhyB6IG51bWVyeWN6bnljaCBixYLEmWTDs3cKcHJ6eSBwb3LDs3dueXdhbml1IGRvIHplcmEuCgpgYGB7ciByZXN1bHRzPUZBTFNFfQoKcGFydF9pZHMgPC0gaGlnaF9kYXRhIHw+IHNlbGVjdChQQVJUX0lEKSB8PiB1bmlxdWUoKQpwYXJ0X2lkcyA8LSBwYXJ0X2lkcyRQQVJUX0lECgojIHBhcmFtcwpkZXRyZW5kID0gRkFMU0UKc2FtcGxpbmcgPSAxMjAKbG93X2ZyZXFfYmFuZCA9IDQKdXBwZXJfZnJlcV9iYW5kID0gMzIKCmZvciAocGlkIGluIHBhcnRfaWRzKSB7CiAgdyA8LSBoaWdoX2RhdGEgfD4KICAgIGZpbHRlcihUcmlhbF90eXBlID09ICJleHBlcmltZW50IiwgUEFSVF9JRCA9PSBwaWQpIHw+CiAgICBncm91cF9ieShDU0kpIHw+CiAgICBzdW1tYXJpc2UobWVhbiA9IG1lYW4oQ29ycikpIHw+CiAgICBhbmFseXplLndhdmVsZXQoCiAgICAgICJtZWFuIiwKICAgICAgbG9lc3Muc3BhbiA9IGRldHJlbmQsCiAgICAgIGR0ID0gMSAvIHNhbXBsaW5nLAogICAgICBsb3dlclBlcmlvZCA9IDEgLyB1cHBlcl9mcmVxX2JhbmQsCiAgICAgIHVwcGVyUGVyaW9kID0gMSAvIGxvd19mcmVxX2JhbmQsCiAgICAgIG1ha2UucHZhbCA9IFRSVUUsCiAgICAgIG4uc2ltID0gMTAKICAgICkKICAKICB3dC5pbWFnZSgKICAgIHcsCiAgICBjb2xvci5rZXkgPSAicXVhbnRpbGUiLAogICAgbWFpbiA9IHBpZCwKICAgIG4ubGV2ZWxzID0gMTAwLAogICAgbGVnZW5kLnBhcmFtcyA9IGxpc3QobGFiID0gIndhdmVsZXQgcG93ZXIgbGV2ZWxzIiwgbWFyID0gNC43KQogICkKICAKICByZWNvbnN0cnVjdCgKICAgIHcsCiAgICBwbG90LndhdmVzID0gVFJVRSwKICAgIGx3ZCA9IGMoMSwgMiksCiAgICBsZWdlbmQuY29vcmRzID0gImJvdHRvbWxlZnQiLAogICAgbWFpbiA9IHBpZAogICkKfQpgYGAK